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ABSTRACT 

We have developed an algorithm for transferring radiation in three- 
dimensional space. The algorithm computes radiation source and sink terms 
using the Fast Fourier Transform (FFT) method, based on a formulation in 
which the integral of any quantity (such as emissivity or opacity) over any 
volume may be written in the classic convolution form. The algorithm is fast 
with the computational time scaling as A^(logA^)^, where N is the number of 
grid points of a simulation box, independent of the number of radiation sources. 
Furthermore, in this formulation one can naturally account for both local 
radiation sources and diffuse background as well as any extra external sources, 
all in a self-consistent fashion. Finally, the algorithm is completely stable and 
robust. 

While the algorithm is generally applicable, we test it on a set of problems 
that encompass a wide range of situations in cosmological applications, 
demonstrating that the algorithm is accurate. These tests show that the 
algorithm produces results that are in excellent agreement with analytic 
expectations in all cases. In particular, radiation flux is guaranteed to propagate 
in the right direction, with the ionization fronts traveling at the correct speed 
with an error no larger than one cell for all the cases tested. The total number 
of photons is conserved in the worst case at ~ 10% level and typically at 1 — 5% 
level over hundreds of time steps. As an added advantage, the accuracy of the 
results depends weakly on the size of the time step, with a typical cosmological 
hydrodynamic time step being sufficient. 
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1. Introduction 

Radiative transfer in three dimensional space is a seven-dimensional problem: three 
spatial dimensions, two angular dimensions, one frequency dimension plus the time 
dimension. As a result, although the basic physics involved is well understood, a direct 
computation of radiative transfer in three dimensional space is prohibitively costly. 

However, radiation field is known to play a very important role in determining the 
ionizational and thermodynamic state of the cosmic gas. Rudimentary treatments of 
radiative transfer from assuming a uniform radiation field (e.g., Cen & Ostriker 1993) to 
using the somewhat improved method with local optical depth approximation (Gnedin & 
Ostriker 1997; Cen & Ostriker 1999) have afforded the breakthrough computations of the 
optically thin regions of the Lya forest with notable successes (Cen et al. 1994; Zhang et 
al. 1995; Miralda-Escude et al. 1996; Hernquist et al. 1996). Nevertheless, fiuctuations in 
the ionizing radiation field are expected to exist even in the optically thin regions, since 
the radiation sources (galaxies and quasars) as well as density fields are known to be 
significantly clustered at all times. Turning away from optical thin regions, one is faced with 
transitional regions of optical depth of order unity and optical thick regions largely shielded 
from external radiation. Moreover, radiation sources themselves most likely reside in dense, 
at least partially shielded regions often with complicated structures and geometries. For all 
these regions a proper treatment of radiative transfer is clearly demanded. The predictive 
power of any theory of galaxy and structure formation, especially with respect to the 
cosmological reionization, the formation and evolution of Lya forest, damped Lya systems 
and galaxies/ stars, is fundamentally limited until one can accurately treat this essential 
process of radiative transfer. 

Significant progress has been made in recent years in developing practical algorithms 
for radiative transfer in cosmological applications. Most of the effort has largely been 
concentrated on two types of approaches: the ray tracing approach (Abel et al. 1999; 
Rozoumov & Scott 1999; Kessel-Deynet & Burkert 2000) and moment equations approach 
(Stone, Mihalas & Norman 1992 for two dimensional space; Norman, Paschos, & Abel 1998; 
Gnedin & Abel 2001). Several Monte Carlo methods have also been explored (Sokasian 
et al. 2001; Ciardi et al. 2001). The primary hmitation of the ray tracing algorithm is its 
intrinsic high cost to follow a large number of radiation sources. The moment equations 
approach condenses the intrinsic difficulty of treating long range radiation propagation out 
as the Eddington tensor. The problem then boils down to evaluating the Eddington tensor 
accurately and efficiently. The recent suggestion by Gnedin & Abel (2001) of computing 
the Eddington tensor by ignoring optical depth is novel and clearly worth being explored 
further. 
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In this paper we develop an entirely new algorithm for transferring radiation in three 
dimensional space. In essence, this algorithm compTites the optical depth between any 
pair of points and the source function for each point using convolution techniques. The 
angular discretization is performed at the receiving site, rather than at the source site 
as in the normal ray tracing scheme, which guarantees coverage of all space regardless 
of the fineness of the angular discretization. The overall scaling of the cost with this 
algorithm is N{\ogNy. Unlike conventional ray casting schemes, the computational cost 
with the present algorithm is independent of the number of sources, suitable for problems 
encountered in cosmological simulations where a large number of direct, ionizing sources 
(e.g., small galaxies) as well as secondary, processed ionizing sources (e.g., scattered 
photons, recombination photons) are present. This paper is organized as follows. In §2 we 
describe the basic algorithm. In §3 we present a battery of tests to show that the present 
algorithm is accurate. Discussion and conclusions are given in §4. 

2. A New Algorithm for Radiative Transfer in Three-Dimensional Space 

In the rest frame of a bundle of photons the equation of radiation transfer for those 
photons (Spitzer 1980) is: 



where s is the path length and Iy{x,n,t) is the specific intensity with I^dududAdt being 
the energy during a time interval dt passing through an area dA about a spatial point 
within a frequency interval dv, within the solid angle duj about t) and ju{x,t) are 

the opacity and emissivity at x at time t. We can integrate equation (1) and write it in an 
integral form for the specific fiux defined as 



dh _ 



(1) 




(2) 



at position x and time t: 




(3) 



where j{x + sn,t) is the emissivity at distance s from position x in the direction n (an 
unit vector) at time t; the two integrals are over path length and solid angle, respectively; 
T,^{x, sn, t) is the optical depth from position x to position x + sn at time t: 
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Tt,{x,sn,t)— / K,y{x, s'n,t)ds' , (4) 
Jo 

where k{x, sn, t) is the opacity at distance s from position x in the direction n at time t. 

The imphcit principal assumption that has been made to derive Equation (3) is that 
the speed of hght is infinity, analogous to the conventional treatment of gravitational 
interactions. Although this assumption may not be necessary, it makes implementation 
of this method especially simple and we adopt it. This assumption is generally excellent 
for cosmological simulations, where light crossing time over a simulation box is indeed 
much smaller than the typical hydrodynamic time step. It should, however, be noted that 
this assumption does not impose any practical restrictions on the propagation speed of 
ionization fronts. 

Equation (3) may be discretized in the following form: 



s n 



where the two sums are over the ray path (assuming to be a straight line) and the solid 
angle, respectively, and the source term S^, is defined as 



Sv{x, sn, t) = jv{xi sn, t)AV{x, sn) (6) 

with AV{x, sn) being the discretization volume element about position x + sn and j^, is the 
mean emissivity over AV{x, sn). The problem now translates to the evaluation of the two 
terms on the right hand side in Equation (5), S,^{x, sn, t) and e"'^"'^^'*"''*). This is the core of 
the overall computation and how it is computed determines the effectiveness and accuracy 
of the algorithm. We propose to evaluate these two terms in a novel way. 

To proceed further we will now choose the spherical coordinate system for the radiation 
field about x, under which we will discretize the angular and radial dimensions in a 
spherically symmetric fashion. This is to say, one can always find a pair of equal volume 
elements (of arbitrary domain shape) one about point x and the other point —x related by 



AV{x,rn) = AV{x, —rn), (7) 

where r is the radial distance from x along a radial direction n or —n. For every vector 
y — rn that lies in AV{x,rn), one can find an opposite vector (about x) —y — r{—n) 
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that lies in AV{x, —rn). This property will be essential in our implementation of the 
computation of the two terms in Equation (5). Let us rewrite the source term S^, in 
Equation (6) as 

Si,{x,rn,t) ^ I ji,{x + rn,t)dV, (8) 

where the integration domain is the volume element AV{x,rn) = r'^drduj{r,n) about 
position X + rn. It is now time to make a simple but algorithmically essential modification 
of Equation (8) by inserting a window function in the integral: 

S„{x,rn,t)^ j j^{x + rn,t)g{r)d^f, (9) 

where g{r) is a window function about the origin x with the following property: 

glr) — 1 for vectors within AV 

— otherwise. (10) 

Note that the integral in Equation (9) is now over the entire space about x. The reader is 
then asked to make the following critical observation: for each window function g'(r) there 
is another window function /(r), which is spherically symmetrically paired with it, i.e., 

/(-r) = 5(r), (11) 
as indicated by Equation (7). Replacing g{r) with f{—r) in Equation (9) gives. 



S^{x,rn,t) = jji,{x + rn,t)f{—r)df. (12) 
Changing integration variable from ftoy = x + r Equation (12) becomes 

S^{x, rn, t) = j ju{y, t)f{x - y)d^y. (13) 

The integration domain is still over the entire space but the window function f{x — y) 
now has a moving origin x. Equation (13) is of the standard convolution form and can be 
efficiently evaluated for all points x simultaneously using the Fast Fourier Transform (FFT) 
techniques, as the Fourier transforms are related by 

s^-kh^ (14) 
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where S_j^, and are Fourier transforms of Si^, jt, and /, respectively. 

The same formulation can be constructed for the optical depth Ti,{x,rn,t). Here we 
instead compute the mean opacity of each volume element 



_ ^ J^yKuix + rn,t)dV 
K^[x,rn,t) = — 

^ I i^u(y,t)f{x -y)d^y 
AV 



(15) 



using the same technique. Then we integrate Ki, radially along n in real space to obtain 
T„{x, rn, t). 



rr 

TJx,rn,t) — / KJx,r'n,t)dr' . (16) 
Jo 

Using the FFT method, Sy{x,rn,t) (Equation 13) and R,^{x,rn,t) (Equation 15) for all 
grid points of a simulation box with respect to all other points in the simulation box can 
be evaluated simultaneously with the total cost scaling as N (log N)'mr'ma, where is the 
total number of grid points, rrir is the number of radial bins along each cone and rUa is the 
number of angular discretization elements. 

The choice of mr and how to discretize the radial direction may be problem dependent. 
It is, however, illuminating to note that both radiation flux and gravitational force 
follow the same inverse square law. It is also noted that the effect due to optical depth 
attenuation tends to relatively suppress the contribution from distant sources. Thus, one 
may conservatively apply some of the techniques and ideas used in gravitational tree codes 
(Barnes & Hut 1986). Specifically, one can more coarsely sample more distant regions to 
discretize the radial direction, the simplest of which would result in oc log N. The choice 
of THa requires some experiments. We note that, in the present algorithm, any coarse level 
of angular discretization guarantees that all sources inside the simulation box are covered. 
This property is in contrast with the conventional ray casting scheme, where a large number 
of rays from a source is necessary in order to intersect every cell in the simulation box at 
least once. The primary operational difference is that in the present algorithm the operation 
is "gathering" flux by the receiving region under consideration, whereas in the ray casting 
method the operation is "broadcasting" flux to each receiving region. Testing shows that 
ma = 16 X 16 = 256 to cover the full solid angle An appears adequate. Combining the three 
scaling factors gives the total number of operations for evaluating Si, or Tj, for all grid points 
at each time step: 
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Total number of operations Ngp — AY-FTN{}.ogNYma, (17) 

where A-p^T is the normal prefactor of the FFT operation. Note that the algorithm 
presented here has the desired property that its computational cost is independent of the 
number of radiation sources. As a matter of fact, each grid point is treated as both a source 
and sink. This property allows us to compute problems where a large number of radiation 
sources are present. Examples include reionization of the universe by sub-galactic galaxies 
and radiation field where diffuse sources such as recombination photons and scattered 
photons contribute significantly. 

We will now implement the algorithm in a three-dimensional simulation box with a 
uniform mesh for hydrodynamic variables and periodic boundary conditions for both hydro 
quantities and radiation field. Two separate coordinate systems and two correspondingly 
separate grids are in use now: the Cartesian coordinate system for the hydrodynamic 
quantities with a uniform Cartesian grid and the spherical coordinate system for the 
radiation field about each hydro grid point with a spherical grid. The angular discretization 
about each hydro grid point is made uniformly across the solid angle 47r with uia elements 
each having a full solid angle of in/rria. The radial discretization is performed using rrir 
bins from r = to Vmax — V^L spaced uniformly in logarithm maximizing the shape of each 
volume element as a cube, where L is the simulation box size. The rrirma window functions 
(i.e., / in Equations 13,15) for each discretized element of the radiation grid in the spherical 
coordinates are computed once initially and their Fourier transforms are stored and used 
for all subsequent time steps. 

So far we have operated in a static coordinate system. In cosmological applications 
one needs to take into account the cosmological effects. We will separate the overall 
radiation field into two parts: 1) radiation fiux due to sources within the simulation box, 
and 2) external radiation fiux. As long as the simulation box is much smaller than the 
Hubble radius, this division allows us to split these two kinds of sources cleanly in the 
simple fashion: radiation fiux from (1) is treated as local, which is not affected by cosmic 
expansion, and radiation flux from (2) is diffuse and subject to the cosmological effects. 
The two kinds of radiation fields are related and computed in a self-consistent way. We 
keep track of radiation that leaves the simulation box and add it to the diffuse background 
(see below) and every point in the simulation box is subject to the diffuse background 
with self-consistently computed attenuation. Thus, the total flux at any point inside the 
simulation box is 
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t) 

{x, t) + t), (18) 

where Fi,^in{x,t) is the contribution from sources inside the simulation box in Equation 
(3) and F^^gxt{^,t) is the contribution from the diffuse background related to the specific 
intensity of the diffuse background at time t, I^^diff{t): 

= h4^ffit)—Ee-^^^'^'''^''''''\ (19) 

where the integral or the sum is over the solid angle about position x; t^{x, Vuppfi, t) is total 
optical depth along the cone about n from position x to a distance of rupp at time t {vupp is 
the upper radius of that cone given the constraint that the simulation box is periodic; Vupp 
is orientation-dependent because the hydro simulation box is not spherical about each grid 
point for the radiation field; the range for r„pp is L < Vupp < Tmax = \/^L). The evolution of 
the specific intensity of the diffuse background, Iy^cLiffif)-i is governed (Cen & Ostriker 1993) 
by 



+ CS„^in{t)+CS^^ext{t)-CKt.,diff{t)I^,diff{t), (20) 

where H{t) is the Hubble constant, c is the speed of light and the mean opacity Rv^diff{t) 
at time t is determined self-consistently by averaging the opacities in the simulation box 
properly taking into account optical shielding effect: 

X 

1 ma 

= — — X^fi;,(f)^e-"''''(^'''"--^'*), (21) 

" X 1=1 

where the outer sum is over all N grid points of the simulation box each having an opacity 
K,^{x) and inner sum is over the solid angles about each point. The term on the right 



dlu4iff{t) 
dt 



^ 0^ •^ii^,diff\t) 



-9- 



hand side of Equation (20) due to sources inside the simulation box is s^^^in, which may be 
computed from emissivity jv{x) of each hydro grid point of the simulation box: 

X 

1 ma 

= ^Ei'^(^)Ee"^'"'^''''"^^"'^ (22) 

" X 1 = 1 

this is the radiation due to internal sources that leaves the simulation box. There may 
be additional sources to the diffuse background. For example, one may wish to take into 
account rare sources such as quasars that are outside the simulation box. This can be 
achieved by simply adding the needed terms as Su,extit)- 

The formulation of Equation (15) deserves more discussion. The integration over the 
domain AV is effectively collapsing the angular dimensions to lump together absorbing 
material into a radial "line" of length Ar for that volume element. This could cause 
an overestimate of opacity for radiation sources that are located inside the same volume 
element or behind the considered volume element but are spatially displaced perpendicular 
to the line of sight with respect to the absorbing material. In a more extreme example, a 
handful of cells of very high optical depth may totally block out radiation downstream, 
if the mean optical depth is evaluated by linearly averaging optical depth spatially. We 
propose to modify the original formulation (Equation 15) as follows. 



^ -K^,{x+sn,t)AL _ JAV^ ^ 

AV 



~ Ay ^^^^ 

and 

R^{x + sn, i) = - — In < e-'^''(*+*"'*)^^ >, (24) 

where AL is the size of a hydro cell. This modified formulation is in some sense like 
averaging the flux instead of the optical depth. This physically motivated modification 
quite successfully circumvents the forementioned possible problems in reahstic situations, as 
subsequent tests will show (see §3). We adopt this particular scheme to present quantitative 
results in the next section (§3). We note that in the case of a uniform opacity Equation 
(15) and Equations (23,24) would yield identical results. 
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Let us briefly summarize the salient features of the proposed algorithm. What we 
have shown is that the algorithm is fast, scaling as A^(logA^)^, where N is the number of 
grid points in a simulation box, independent of the number of radiation sources. Given 
that computation does not involve any differentiation and all computed integrals are 
guaranteed to converge in all situations, the algorithm is completely stable and robust. We 
emphasize that the proposed algorithm, by design, guarantees that flux propagates in the 
right direction in any situation. The tests presented in the next section will show that the 
amplitude of flux is also computed with very small errors, indicating that the algorithm is 
also accurate. 



3. Tests of the Algorithm 

We will present a suite of tests, which are intended to have significant bearings on 
cosmological applications. We will show that the method performs very well in all tested 
situations. Clearly, fine tuning of the scheme is worth exploring and is deferred for future 
work. 

For the tests shown we use a 32'^ Cartesian grid with periodic boundary conditions. 
We use rUa = 256 uniformly spaced angular bin and = 18 logarithmically spaced radial 
bins for all cases. For the tests shown we do not include cosmological affects; i.e., H{t) is 
set to zero. But the accuracy of the algorithm should not be affected by this simplification. 



3.1. Static Spherical Optical Depth Distribution 

It is guaranteed that the proposed formalism will give the correct flux distribution in 
the case of no optical attenuation. The simplest non-trivial test then is a source embedded 
in a static, uniform, non-zero opacity distribution. The flux distribution about a single 
source is then 



= i^^""' (25) 

where L is the luminosity of the source; r is source-centric radius; R is the opacity. Figure 
ID shows the flux distribution in the x-y plane with z = 16, where the luminous source is 
located at cell (16, 16, 16) and the opacity per cell is R = 0.5. The solid contours nearly 
overlap with the analytic expectations (dashed contours) with some slight deformations 
from sphericity simply due to the nature of the Cartesian grid. Clearly, the computed flux 
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is accurate. 

Since the equation for the flux is hnear with respect to the source luminosity, the 
method would obviously yield results of the same accuracy for any number of distributed 
sources. 



3.2. Static Elliptical Optical Depth Distribution 

The next step is to test a non-spherical static opacity distribution. An axisymmetric 
elliptical distribution of opacity of the following form is tested: 

n{e) = (26) 

'1- (1-62)C0S2^ 



where 6 is the angle with respect to the positive z direction, and c and b are the 
normalization and the ellipticity parameter of the distribution, respectively; k = c at 
cos 6' = 1 and k = cb at cos 6' = 0. 

Figure shows the computed flux distribution (solid contours) in the x-z plane with 
y = 16 compared to the analytic result (dashed contours). We use c = 0.2 and b = 0.2 (see 
Equation 26) for this illustration. We see that the results are quite satisfactory given the 
Cartesian grid. The agreement with the analytic result is very good except near cos^ = 1 
where the opacity is the highest and a maximum error of 2 cells is made there. This is 
caused by a slight underestimate of the optical depth there, due to the averaging scheme 
of optical depth as indicated by Equations (23,24). We have also tested similar cases by 
varying either the opacity c or ellipticity parameter b (see Equation 26) and find that the 
results have similar accuracies. 



3.3. Ionization of a Uniform Neutral Medium 

We now switch to a higher gear to test the algorithm in a situation where the 
distribution of opacity changes self-consistently with time due to ionization. For this and 
subsequent tests we monitor photon number conservation by computing 

A^e(^) + Ndiff(t) , , 

Vit) = (27) 
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Fig. 1. — Distribution of flux in the x-y plane (with z = 16) for a source sitting at 
cell (16,16,16) embedded in a static, uniform opacity distribution. The contour levels are 
logarithmically spaced with an increment of 1 dex per contour level. The dashed contours 
are the analytic results (Equation 25) and solid contours are obtained with the present 
algorithm. The units for the contours levels are arbitrary. 
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Fig. 2. — Distribution of flux in the x-z plane (with y = 16) for a source sitting at cell 
(16,16,16) embedded in a static, elliptical opacity distribution (Equation 26). The contour 
levels are logarithmically spaced with an increment of 0.5 dex per contour level. The dashed 
contours are the analytic results and solid contours are obtained with the present algorithm. 
The units for the contours levels are arbitrary. 
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where Ne{t) is the number of free electrons created by time t (assuming 100% hydrogen) 
in the simulation box and NfUfj is the number of photons in the diffuse background in 
the simulation box volume at time t and Nph^tot is the total number of photons emitted 
by time t by sources in the simulation box. If photon number conservation is strictly 
observed, t] would be unity. In all the following tests recombination time is set to infinity 
to make the problem simple to track. But the accuracy of the method is unchanged by this 
simplification. 

We start with the case of a uniform neutral medium, consisting entirely of hydrogen 
atoms. An ionizing source of ionizing luminosity Nph photon/sec is placed at cell (16, 16, 16). 
The ionizing photons from the source ionize the neutral medium outward and the radius 
of the ionization front, the surface that separates the interior ionized medium from the 
exterior neutral medium, evolves with time approximately as: 



r{t) = 13.7(^)1/3^ ;L.^)-i/3(^^)i/3kpc, (28) 

where n is the number density of the neutral hydrogen and t is the elapsed time. The reason 
that this formula (26) is only approximate is that some photons may penetrate further out 
ahead of the ionization front, if the optical depth per cell is not sufficiently high, and then 
the "front" is no longer well defined; the formula would be highly accurate for a high opacity 
distribution. The optical depth per cell is = (ThuAx = 15.4(?2/10~3cm~3)(Ax/5kpc), 
where the hydrogen photo-ionization cross section an = 10~^^ cm^ is used throughout and 
Ax is the cell size. Thus, for the present test, the optical depth per cell is very large and 
formula (26) should be very accurate. 

Figure ^ shows the contours of the neutral hydrogen fraction in the x-y plane with 
2; = 16 at four epochs. We use Nph = 10^^ photon/sec, n = 10~^ cm~^ and cell size of 
Ax = 5 kpc. The adopted gas density is about 15 times the mean gas density at redshift 
z = 6. We see that the agreement between the computed results and analytic expectations 
is excellent at all times. The discrepancy on the radius of the ionization front is evidently no 
larger than one cell. Given the approximations that have been made in the implementation, 
especially with respect to the computation of mean opacity (Equations 21,22), it is unclear 
how well the algorithm conserves photons. Figure ^ shows the degree of conservation of 
photons as a function of time. It is seen that, except for the first 8 time steps, the number 
of photons is conserved at better than 4% with the average at about 2 — 3%. Timesteps of 
10'^ yrs are used but the results are insensitive to the timestep; using a timestep 3 times 
larger yield very similar results except for the upper left corner panel at the beginning of the 
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Fig. 3. — Distribution of neutral hydrogen fraction in the x-y plane (with z = 16) for a source 
of luminosity 10^^ photon/sec sitting at cell (16,16,16) in a uniform density n = Ix 10^^ cm^'^ 
at four epochs (5 x 10^, 1 x 10®, 4 x 10®, 1 x 10^) yrs. The solid contours indicate the neutral 
hydrogen fraction of 0.3, 0.6, 0.9 inside out computed with the present algorithm. The dashed 
contours are the radii from analytic calculations (Equation 28). 
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Fig. 4. — The error on photon number conservation (see Equation 27) as a function of time 
for the case of an ionizing galaxy surrounded by a uniform density distribution (see Figure 
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simulation. Note that the age of the universe at z = 6 is about 5 x 10^ yrs. The indicated 
convergence of results at this relatively large size of timesteps is clearly advantageous. 



3.4. Ionization of an r ^ Neutral Medium 

Next, we examine a case of a source sitting in a spherical density distribution, 
n = nc{r/rc)~^ cm~^; the central cell at r = is set to a density equal to nc{Ax/rc)~'^, 
where Ax is the cell size. This density profile is close to that in the central regions of dark 
matter halos (Navarro, Frenk, & White 1997) and is of significant cosmological relevance. 
The expected evolution of the radius of the ionization front is approximately: 



r{t) 



^ 35.2( J ( + 25 (— )2/3 - — ^ kpc,(29) 

\ ^105isec-i^^l0-2cm-3^ ^lOVs 5kpc^ VA'k' 2T^)\h\iY>cj ^ ' 



where the second (correction) term inside the square root is due to the non-singular 
core of the modified density distribution. The optical depth per cell as a function of 
source-centric radius is r(r) = 154(ri/10~^cm~^)(Ax/5kpc)(r/rc)~"'^. Therefore, for the 
adopted values of Tc, ric and Ax (see below), the formula will be very accurate over the 
range of radii examined and the thickness of the ionization front should be thin and close 
to one cell. 

Figure ^ shows the contours of the neutral hydrogen fraction in the x-y plane with 
z = 16 at four epochs. We use Nph = 10^^ photon/sec, ric = 1.5 x 10~^ cm~^, Tc = 5 kpc and 
cell size of Ax = 5 kpc. Note that the adopted ric is approximately 200 times the mean gas 
density at z = 6. Timestep used is 10^ yrs but the results are insensitive to the timestep. 
From Figure ^ it is seen that the agreement between the computed results and analytic 
expectations is excellent at all times and the discrepancy on the radius of the ionization 
front is no larger than one cell. Figure |^ presents the degree of conservation of photons as a 
function of time, showing that the number of photons is conserved at better than 10% with 
the average at about 2 — 3% in this case. 



3.5. Ionization of an r ^ Neutral Medium 

We turn to a steeper density distribution, n = ric^r/rc)'"^ cm~^; the central cell at 
r = is set to a density equal to nc(Ax/rc)~^, where Ax is the cell size. The density 
profile is close to that of halos just interior of the virial radius (e.g., Tyson, Kochanski, 
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Fig. 5. — Distribution of neutral hydrogen fraction in the x-y plane (with z = 16) for 
a source sitting at cell (16,16,16) in a density distribution of n(r) oc at four epochs 
(5 X 10'', 2 X 10*^,5 X 10^,1.5 x 10^) yrs. The dashed contours are the analytic results 
(Equation 29) and solid contours are obtained with the present algorithm indicating the 
neutral hydrogen fractions of 0.3, 0.6, 0.9 inside out. 
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Fig. 6. — The error on photon number conservation as a function of time for the case of an 
ionizing galaxy surrounded by an density distribution (see Figure 
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& dell' Antonio 1998). The expected evolution of the radius of the ionization front is 
approximately: 

where the second term on the right hand side is due to the non-singular core imposed on 
the density profile. 

Figure |^ shows the contours of the neutral hydrogen fraction in the x-y plane with 
z = 16. We use Nph = 10^^ photon/sec, ric = 1.5 x 10~^ cm~^, = 5 kpc and cell size of 
Ax = 5 kpc at four epochs. Timestep used is 10'' yrs but like in other cases the results 
are insensitive to the timestep. The computed results and analytic expectations agree very 
well at all times. For the adopted values of Uc, Tc and Ax, the optical depth per cell as a 
function of source-centric radius is r(r) = 154(?7,/10~^cm~^)(Ax/5kpc)(r/rc)~^ shown in 
Figure we have r(r) = 3.6 for r = 8 cells. Thus, the computed shell of the ionization front 
will be somewhat broadened at late times due to a deeper penetration of ionizing photons 
outward of the ionization front, clearly evident in Figure |^ Figure ^ shows that the number 
of photons is conserved at better than 6% with the average at about 1 — 4% for this case. 



3.6. Ionization of an Elliptical r ^ Neutral Medium 

In realistic situations a source may sit inside a gas cloud that may not be spherical. 
We examine an elliptical isothermal distribution with the following density distribution 



n{r,9) = nc{r/rc) ^ , (31) 

a-(l-62)cos2^ 



where 6 is the angle with respect to the positive z direction, and b indicates the ellipticity 
of the distribution. The central cell at r = is designed not to be singular but has a density 
equal to nc(Ax/rc)~^, where Ax is the cell size. The expected, orientation-dependent 
evolution of the radius of the ionization front follows Equation (30). 

Figure |l^ shows the contours of the neutral hydrogen fraction in the x-z plane with 
?/ = 16 at four epochs. We use Nph = 10^^ photon/sec, Uc = 1.5 x 10~^ cm~^, Tc = 5 kpc 
and cell size of Ax = 5 kpc. Timestep used is 1.3 x 10^ yrs but the results are insensitive to 
the timestep. It is seen that the algorithm handles the anisotropic density distribution quite 
well and the agreement between the computed results and analytic expectations is very 
good at all times, with the computed ionization "fronts" tracing out the analytic results 
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Fig. 7. — Distribution of neutral hydrogen fraction in the x-y plane (with z = 16) for 
a source sitting at cell (16,16,16) in a density distribution of n(r) oc at four epochs 
(1.2 X 10'', 4 X 10'', 8 X 10'^, 2 x 10*^) yrs. The dashed contours are the analytic results 
(Equation 30) and solid contours are obtained with the present algorithm indicating the 
neutral hydrogen fractions of 0.3, 0.6, 0.9 inside out. 
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Fig. 8. — The original distribution of optical depth as a function of radius with an 
density distribution. We see that at r > 6 cells a substantial amount of photons starts to 
penetrate outward ahead of the ionization front causing the slight thickening of the front, 
see in Figure |^. 
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Fig. 9. — The error on photon number conservation as a function of time for the case of an 
ionizing galaxy surrounded by a density distribution (see Figure |^). 
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Fig. 10. — Distribution of neutral hydrogen fraction in the x-z plane (with y = 16) 
for a source sitting at cell (16,16,16) in an elliptical isothermal density distribution of 
n{r) oc r~'^f{6) [see Equation 31 for the elliptical function of f{0)] at four epochs, 
(1.3 X 10'^, 4 X 10'', 8 X 10'', 1.2 x 10®) yrs. The dashed contours are the analytic results 
and solid contours are obtained with the present algorithm, indicating the neutral hydrogen 
fractions of 0.3, 0.6, 0.9 inside out. 
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Fig. 11. — The error on photon number conservation as a function of time for the case of 
an ionizing galaxy surrounded an elhptical isothermal density distribution of n(r) oc r~'^f{9) 
(see Figure |l^) . 
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quite nicely with an error no larger than about one cell. The thickening of the computed 
"fronts" is real due to penetration of photons. Figure |TT] shows the degree of conservation 
of photons. It is seen that the number of photons is conserved at better than 6% with the 
average at about 1 — 4%. 



3.7. Ionization of a Sharply Divided Anisotropic Neutral Medium 

So far we have examined several different cases with smooth density distributions. It is 
prudent to test the algorithm in a more challenging situation. We pick a case where there 
are two different density fields separated by a sharp boundary. An axisymmetric density 
field is prescribed as 

n(e) = lO-^cm-^ for - < e < — 

4 - - 4 

= 10~'^cm~'^ otherwise. (32) 

The source sits at cell (16, 16, 16). The expected evolution of the radius of the ionization 
front follows Equation (28) for each of the two domains with different speeds of the 
ionization fronts. 

Figure shows the contours of the neutral hydrogen fraction in the x-z plane with 
y = 16 at four epochs. We use Nph = 10^^ photon/sec, Tc = 5 kpc and cell size of 
Ax = 5 kpc. Timestep used is 6.7 x 10^ yrs but the results are insensitive to the timestep. 
The agreement between the computed results and the analytic results are surprisingly 
good with the radius of the ionization front accurate to about one cell. At the boundaries 
separating the two regions there is a slight over-shielding on the low density side affected by 
the high optical depth, high density side. The width of this "affected" region is limited by 
the size of the solid angle of each discretized angular element, which in the present case is 
about 12^ square degrees. Figure |13| shows the degree of conservation of photons, indicating 
that the number of photons is conserved at better than 10% with the average at about 
1 - 5%. 



3.8. Double Sources in a Uniform Density 

Up to now we have tested cases with a single radiation source. It is important that 
the algorithm is capable of handling interactions of multiple radiation sources through 
the process of percolation, a process that is thought to occur during the cosmological 
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Fig. 12. — Distribution of neutral hydrogen fraction in the x-z plane (with y = 16) for a 
source sitting at cell (16,16,16) embedded in a sharply divided anisotropic density field (see 
Equation 32) at four epochs (2 x 10^, 2 x 10'', 4 x 10'', 7.3 x 10'^) yrs. The dashed contours are 
the analytic results and solid contours are obtained with the present algorithm, indicating 
the neutral hydrogen fractions of 0.3, 0.6, 0.9 inside out. 
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Fig. 13. — The error on photon number conservation as a function of time for the case of an 
ionizing galaxy surrounded by a sharply divided anisotropic density distribution (see Figure 

m 
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Fig. 14. — Distribution of ionization fraction in the x-y plane (with z = 16) for two 
static sources sitting at cell (23,16,16) and cell (13,16,16) with luminosities of Nph = 
2 X 10^° photon/sec and Nph = 10^^ photon/sec, respectively, with a separation of 50 kpc, 
at four epochs (5 x 10^, 1 x 10*^, 2 x 10®, 7 x 10®) yrs. The dashed contours are the analytic 
results and solid contours are obtained with the present algorithm, indicating the neutral 
hydrogen fractions of 0.3, 0.6, 0.9 inside out. 
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Fig. 15. — The error on photon conservation as a function of time for the case of two ionizing 
galaxies with luminosities of 10^^ sec~^ and 2 x 10^° sec~^, respectively, separated by 50 kpc 
surrounded by a uniform density distribution with n = 1 x 10^'^ cm~^. 
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reionization. We will demonstrate this by studying the situation where there are two sources 
of unequal luminosities. The two sources with luminosities Nph = 10^^ photon/sec, and 
2 X 10^° photon/sec, respectively, sit in a uniform neutral medium of density n = 10~^ cm~'^ 
initially. A cell size of 5 kpc is used for the simulation box. The two sources are located at 
cell (23,16,16) and cell (13,16,16) with a separation of 10 cells (i.e., 50 kpc). Timestep used 
is 10^ yrs but the results are insensitive to the size of the timestep. 

Prior to the overlap of the two HII regions produced by the two sources separately, 
the evolution of each HII region is separate and the evolution of the radii of the ionization 
fronts follows Equation (28). Subsequent to the overlap their combined HII region will 
continue to expand but an analytic solution is not readily available. Figure |n| shows the 
contours of the neutral hydrogen fraction in the x-y plane with z = 16 at four epochs. It 
is clear that the algorithm nicely treats the two separate regions without any interference 
before the overlap, with the two separate ionization fronts traveling at the correct speeds. 
For the post-overlap era we only show the analytic solution as if overlap has not occurred 
for the sole purpose of guiding the reader's eye. The evolution of the ionization front in 
the post-overlap era is more complicated but the computed results appear to be quite 
reasonable, with the overlapped region continuing to expand and becoming rounder with 
time. In this case, one is more keen to check if the total number of photons is conserved. 



Figure 15 shows the degree of conservation of photons. It is seen the number of photons is 



conserved at better than 9% with the average at about 1 — 2%. 



3.9. Quadruple Sources in a Uniform Density 

Let us make the situation a bit more interesting. We will add two more sources to the 
case tested in §3.8. The two additional sources have luminosities Nph = 3 x 10^° photon/sec, 
and Nph = 5 x 10^° photon/sec, sitting at cells (20,8,16) and (8,24,16) but turned on with 
lags of 1 X 10^ yrs and 2 x 10^ yrs (relative to the turn-on time of the initial two sources), 
respectively. 

Figure |1^ shows the contours of the neutral hydrogen fraction in the x-y plane with 
2; = 16 at four epochs. For the post-overlap era we only show the analytic solution as if 
overlap has not occurred for the sake of illustration. We see that the algorithm follows the 
separate HII regions, as expected. Continuous formation of galaxies in time can evidently 
be properly followed. Subsequent mergers of HII regions occur as expected. Figure |l^ 
shows the degree of conservation of photons. It is seen the number of photons is conserved 
at better than 9% with the average at about 1 — 3%. The evolution of the ionization front 
in the post-overlap era is very complicated but the computed results look reasonable, with 
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the overlapped region continuing to expand in a fashion that is expected. The fact that 
photon number is well conserved and flux is designed to travel in the right direction ensures 
that the computed results should be accurate. 



3.10. Shadowing by a Dense Gas Cloud 

We will now step back to a case with a single source. But we will include in the 
uniform density field a dense gas cloud, which is designed to be optically thick to cast a 
shadow in the direction of radiation propagation. It is expected that shadowing should 
be common in the real universe since galaxies and other dense clouds (such as damped 
Lyman alpha systems) are (at least partially) opaque to radiation. A spherical gas cloud 
of radius 20 kpc with a density distribution of 1 x 10"^(r/20kpc)~^ cm~^ is centered at 
cell (24, 16, 16). The rest of the simulation box has a uniform density of 1 x 10~^ cm~^. 
The ionizing source is located at cell (16,16,16) with a luminosity of 10^^ sec^^. A cell size 
of Ax = 5 kpc is used for the simulation. Figure |1^ shows the neutral hydrogen fraction 
contours at four epochs. The analytic results are only valid before the ionization front 
reaches the halo; at subsequent times regions behind the halo within the region delimited 
by the two dotted lines should remain neutral. We indeed see the shadow cast by the dense 
gas cloud. The slightly overshadowing near the edges of the shadowed region is due to the 



limited resolution of angular discretization, also seen in Figure |12 



The relatively complicated situation here does not guarantee that photon number 
conservation will be obeyed. Figure |19| shows the degree of conservation of photons. Quite 
reassuringly, we see that the number of photons is conserved at better than 10% with the 
average at about 1 — 2%. 



3.11. Outside-In Ionization of a Spherical Isothermal Gas Cloud 

Let us next examine another class of ionization processes, where external diffuse 
radiation ionizes isolated overdense regions. This type of ionization process is common in 
cosmological situations where low density regions become ionized first. In this case we check 
the photon number conservation in a slightly different way as: 

Vit) = — -, (33) 

" incoming \ J 

where Ne{t) is the number of free electrons created by time t in the simulation box and 
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Fig. 16. — Distribution of neutral hydrogen fraction in the x-y plane (with z = 16) at 
four epochs (5 x 10''', 1 x 10^, 2 x 10^, 7 x 10^) yrs. for four sources sitting at cells (23,16,16), 
(13,16,16), (20,8,16) and (8,24,16), with luminosities of iVp;, = (2x10^°, 105\ 3x10^°, 5x10^°) 
photon/sec, with turn-on times at t = (0,0, 1 x 10^,2 x 10^) yrs, respectively. The dashed 
contours are the analytic results (only valid times before the HII regions overlap) and solid 
contours are obtained with the present algorithm, indicating the neutral hydrogen fractions 
of 0.3, 0.6, 0.9 inside out. 
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Fig. 17. — The error on photon number conservation as a function of time for the case of four 
ionizing galaxies with luminosities of Nph = (2 x 10^°, 10^^,3 x 10^°, 5 x 10^°) photon/sec, 
with turn-on times at t = (0, 0, 1 x 10*, 2 x 10*) yrs, respectively, surrounded by a uniform 
density distribution with n = 1 x 10~^ cm~^ (see Figure |1^). 
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Fig. 18. — Distribution of neutral hydrogen fraction in the x-y plane (with z = 16) for a 
source sitting at cell (16,16,16) of luminosity 10^^ sec~^ surrounded by a uniform density 
distribution with n = 1 x 10^^ cm~^ at four epochs (4 x 10'', 1 x 10®, 5 x 10®, 1 x 10^) yrs. 
In addition, there is a spherical gas cloud of radius 20 kpc centered at (24, 16, 16) with a 
density distribution of 0.1(r/20kpc)^^. The small dashed circle on the right in each panel 
indicates the size of the halo. The (larger) dashed contours (on the left) are the analytic 
results and solid contours are obtained with the present algorithm, indicating the neutral 
hydrogen fractions of 0.3, 0.6, 0.9 inside out. The analytic results are only valid before the 
ionization front reaches the halo; at subsequent times regions behind the halo delimited by 
the two dotted lines should remain neutral and the (large) circle on the left serves only to 
indicate the position of the ionization front if there were no obscuration. 
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Fig. 19. — The error on photon number conservation as a function of time for the case of one 
ionizing galaxy with luminosity of 10^^ sec~^ surrounded by a uniform density distribution 
with n = 1 X 10"'^ cm~^. In addition, there is a spherical gas cloud of radius 20 kpc centered 
at (24, 16, 16) with a density distribution of 1 x 10"^(r/20kpc)~^ 
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Nincoming IS the number of incoming photons from the diffuse background that are absorbed 
in the simulation box volume by time t. If photon number conservation is strictly observed, 
7] would be unity. Recombination time is set to infinity. 

The first simple example is the ionization of a spherical isothermal gas cloud of density 
n = nc{r/r^~'^ cm~^ with an outer cutoff radius r^; the central cell at r = has a density 
equal to ?T,c(Ax/rc)"^ to avoid singularity, where Ax is the cell size. The expected evolution 
of the radius of the ionization front, which now travels inward towards the center of the gas 
cloud, is approximately: 

2 

r{t) = -, (34) 

where F is the diffuse flux assumed to be constant with time. Note that the factor 0.5 in 
front of flux F in Equation (34) is due to the fact that any element on the surface is subject 
to radiation only at half the total solid angle. 

Figure shows the evolution of the ionization front in the x-y plane with z = 16. We 
use Tc = 5 kpc, Uc = 1.5 x 10~^, Ax = 1 kpc, b = 0.5 and F = 1.9 x 10^ cm~^sec~^. We see 
that the analytic evolution of the ionization front is well tracked by the computed results 
with error on the radius no larger than one cell. The spiky surface at the earliest time 
shown reflects the initial condition laid out on a Cartesian grid. Figure BT] shows the photon 



number conservation as a function of time and again indicates that the method conserves 
total number of photons very well at < 1%. The high degree of photon conservation in case 
is attributable to the fact that there is no obscuration for the cells on the ionization surface 
that receive ionizing photons. 



3.12. Outside-In Ionization of an Elliptical Isothermal Gas Cloud 

We next consider ionization of an isolated elliptical isothermal gas cloud of the following 
density profile by a diffuse ionizing background 



n{r,6) = nc{r/rc) ^ , ; (35) 

- (l-62)cos2^ 



the central cell at r = has a density equal to ?T,c(Ax/rc)~^ to avoid singularity, where Ax 
is the cell size. The expected evolution of the radius of the ionization front follows Equation 
(34) with an angle-dependent Uc indicated by Equation (35). 
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Fig. 20. — Distribution of neutral hydrogen fraction in the x-y plane (with z = 16) for 
an isothermal gas cloud sitting at cell (16,16,16), exposed to an external isotropic diffuse 
radiation background, at four epochs (2 x 10^, 2 x 10^, 6 x 10^, 2 x 10®) yrs. The dashed 
contours are the analytic results (Equation 34) and solid contours are obtained with the 
present algorithm, indicating the neutral hydrogen fractions of 0.3, 0.6, 0.9 outside in. 
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Fig. 21. — The error on photon number conservation as a function of time for the case of 
a halo of density ric = 1.5 x 10^^(r/5kpc)~^ cm~^, subject to a diffuse background flux of 
19 X 10^ cm~^sec"^ 
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Fig. 22. — Distribution of neutral hydrogen fraction in the x-z plane with y = 16 for an 
elliptical isothermal gas cloud (see Equation 35) sitting at cell (16,16,16), exposed to an 
external isotropic diffuse radiation background, at four epochs (2 x 10^, 2 x 10^, 6 x 10'', 2 x 
10^) yrs. The dashed contours are the analytic results and solid contours are obtained with 
the present algorithm, indicating the neutral hydrogen fractions of 0.3, 0.6, 0.9 outside in. 
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Fig. 23. — The error on photon number conservation as a function of time for the case of 
an eUiptical isothermal gas cloud (see Equation 35) sitting at cell (16,16,16), exposed to an 
external isotropic diffuse radiation background. 
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Figure |2^ shows the evolution of the ionization front in the x-z plane with y = 16. We 
use = 5 kpc, ric = 1.5 x 10~^, Ax = 1 kpc and F = 1.9 x 10^ cm~^sec~^. We see that the 
analytic evolution of the ionization front is well traced by the computed results with error 
on the radius no larger than one cell. Figure ^ shows the photon number conservation as 
a function of time and again indicates that the method conserves total number of photons 
very well at < 1%. 



3.13. Radiation Propagation In a Realistic Cosmological Density Field 

Finally, we turn to the radiation propagation in a realistic density field produced by 
cosmological simulations. In this case, the situation is much more complicated and no 
analytic solution is possible. But many of the cases tested above have significant bearings 
on this realistic case. We take a subbox from one of our latest simulations (Cen et al. 2001) 
at redshift z = 6. The subbox has a 32^ grid and a cell size of 4.65 kpc proper. 

We turn on by hand a galaxy of ionizing luminosity of Nph = 10^^ sec~^ at cell 
(11,22,21) corresponding to the densest cell in the box, as indicated by the solid dot in 
Figures (24,25,26). We assume that the gas is composed entirely of hydrogen and at the 
starting time hydrogen is entirely neutral. For simplicity we have set the recombination 
time to infinity. It is noted that the galaxy sits in a generic filament, which is embedded in 
a complex structure containing voids and other interesting structures. Since the subbox is 
a small region of a large, periodic box, it is not periodic itself. But we treat it as if it were 
periodic. 



Figure 53 and Figure 53 show the distribution of neutral hydrogen at time 



t = {8 X 10'', 4 X 10^)yrs (middle columns), respectively, compared to that at the starting 
time {t = 0; left column). The most notable feature is that the radiation propagation is 
highly anisotropic. The ionizing photons quickly clear out widening tunnels towards low 
density (void) regions in the direction perpendicular to the filament, while radiation front 
travels at much slower speeds in other directions. Although this is not at all surprising, 
the figures do remind us of the complex situations encountered in cosmological structures 
and re-iterate the need for a proper treatment of radiative transfer. We also use a high 
resolution ray tracing code to obtain numerically the approximate "true" distribution, 
shown as the right columns. Note that in the ray tracing code one just counts photons and 
compares to the number of hydrogen atoms along each ray cone; photon penetration ahead 
of the ionization front is not taken into account and ionization front is taken to be abruptly 
sharp. Nevertheless, the results obtained with the ray tracing method should be fairly close 
to the 'truth". We see that the agreement between the result obtained with the present 




Fig. 24. — The left column shows projected gas density fields along three faces of the box 
(three rows) for a three-dimensional box with 32'^ grid points and a cell size of 4.65 kpc proper 
at the initial time. The solid dot indicates the position of the putative galaxy of luminosity 
10^^ sec~^. The contours are: the first contour is at the mean density and subsequent ones 
increase by 0.2 dex per contour. The middle column shows projected neutral hydrogen 
density fields along three faces of the box at t = 8 x lO^yrs, computed with the present 
method. The right column shows projected neutral hydrogen density fields along three faces 
of the box at t = 8 X 10''yrs, computed using a high resolution ray tracing code (with 38416 
rays from the source). 
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Fig. 25. — The left column shows projected gas density fields along three faces of the box 
(three rows) for a three-dimensional box with 32'^ grid points and a cell size of 4.65 kpc proper 
at the initial time. The solid dot indicates the position of the putative galaxy of luminosity 
10^^ sec~^. The contours are: the first contour is at the mean density and subsequent ones 
increase by 0.2 dex per contour. The middle column shows projected neutral hydrogen 
density fields along three faces of the box at t = 4 x lO^yrs, computed with the present 
method. The right column shows projected neutral hydrogen density fields along three faces 
of the box at t = 4 X lO'^yrs, computed using a high resolution ray tracing code (with 38416 
rays from the source). 
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Fig. 26. — Volume-weighted flux distributions projected along three faces of the box at the 
two times, t = (8 x 10^, 4 x 10*)yrs, shown in Figure and Figure ^ respectively. The first 
contour around the solid dot (the galaxy) has a flux of 10^'^cm~^sec~^ and the decrement 
for the successive contours is 0.25 dex. 
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Fig. 27. — The error on photon number conservation as a function of time for the case of 
ionization of a reahstic density field produced in cosmological simulations at z = 6. 
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method and that with the ray tracing method is excellent. 

Figure shows the projected flux distributions (volume-weighted) at the two epochs 
respectively. We see that flux distributions are quite complex but consistent with the 
neutral hydrogen distributions. It is seen at the earlier epoch, when substantial blocking or 
shadowing exists, the contours of the flux distribution tend to orient perpendicular to the 
contours of the neutral hydrogen density. Note that in the present calculations we have set 
the recombination time to inflnity; in a more involved calculation including recombination 
it is expected that fluctuations in the flux distributions would be enhanced. 

Finally, Figure shows how the total number of photons is conserved. We are once 
again seeing satisfactory conservation of photons in this realistic simulation at about 1 — 10% 
level with a time averaged accuracy of 3 — 5%. This level of photon number conservation 
is quite adequate for cosmological simulations, since other uncertainties, including galaxy 
luminosity, radiation escape fraction from the galaxies, star formation efficiency, etc., are 
usually much larger. Besides, we expect that we may be able to improve the accuracy of 
the algorithm by more detailed experimenting and careful fine tuning, if necessary. 

The fact that photon number is well conserved and flux is designed under the algorithm 
to propagate in the right directions guarantee that the results should be reliable, as verifled 
by the comparison with the results from a high resolution ray tracing method. 

3.14. Summary on the Tests 

We have performed a variety of tests on the proposed algorithm for transferring 
radiation in three dimensional space. The tests are designed to be relevant to cosmological 
applications. It is worthwhile to re-emphasize that the proposed algorithm, by design, 
guarantees that flux propagates in the right direction in any situation. Thus, the tests 
should mainly demonstrate how accurately the amplitude of the flux is computed. 

To summarize the results we flnd, the flux is computed very accurately, resulting in 
ionization fronts that travel at the correct speed with error no larger than one cell in all 
cases. Photon number is conserved with a maximum error of about 10% and an average 
error of 1 — 5% over hundreds of time steps for all the cases tested, with = 256 angular 
elements and = 18 logarithmically spaced radial bins on a 32^ grid. We note that 
in the limit of inflnite values for and m^, the present method solves the equation of 
the radiative transfer exactly and the approximation would represent the truth. In other 
words, the tests have basically shown that with reasonable number of angular and radial 
discretization elements the proposed method gives quite accurate results. 
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Finally, we find that the accuracy of the algorithm does not depend sensitively on the 
time step used. A time step of size comparable to the typical time step in cosmological 
hydrodynamic simulations appears sufficient. This feature is desirable with regard to the 
overall computational cost. 

4. Discussion and Conclusions 

In light of the recent observations of high redshift quasars (e.g., Fan et al. 2001), 
indicating that the cosmological reionization may be just ending at redshift z 6 
(Becker et al. 2001; Barkana 2001; Cen & McDonald 2001), a new and exciting frontier 
is forming. Inhomogeneous reionization of the imiverse is now within the observational 
reach and could potentially provide a great tool to test both cosmological models and 
astrophysics of galaxy formation (e.g., Barkana & Loeb 2001). It has thus become 
mandatory to perform cosmological hydrodynamic simulations that tackle this problem 
with an adequate treatment of radiative transfer. Furthermore, even in the post-reionization 
era, an initial inhomogeneous reionization process coupled with the significant clustering 
of radiation sources (galaxies or quasars) as well as the fiuctuating density field hkely 
produces fluctuations both in the thermodynamic properties of the cosmic gas and in 
the ionizing radiation field. Such fluctuations may substantially alter quantitatively (or 
even qualitatively) the picture of the optically thin regions of the Lja forest on a variety 
of scales, which would limit our ability of accurately extracting important cosmological 
information from observations of the Lya forest. Lastly (but not the least), since almost all 
neutral gas resides and star formation occurs only in dense, optically thick regions, many 
questions pertaining neutral gas and star formation in the universe can be answered with 
confldence only when radiative transfer is properly included. 

We have developed a fast, accurate and robust algorithm for radiative transfer in 
three-dimensional space. The core of the algorithm is based on a formulation that the 
summation of any quantity (such as emissivity or opacity) over any volume can be written 
in the standard convolution form, which can be computed efficiently using the Fast Fourier 
Transform techniques. We will name this method "FRT" (Fast Radiative Transfer) method. 
The overall computational time with this algorithm scales as A^(logA^)^, where N is the 
number of grid points in a simulation box. Therefore, this is a fast algorithm. We stress that 
the computational cost with this method is independent of the number of radiation sources, 
providing a natural match to cosmological simulations where a large number of sources are 
present. Local sources and diffuse background are naturally split and the evolution of each 
of the two components as well as interactions between them are handled self-consistently. 
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The devised integral form of the algorithm guarantees that the method is completely stable 
and robust. 

The algorithm is tested on a wide range of analytically tractable problems that have 
significant bearings on cosmological applications. We find that the algorithm performs very 
well in all cases, accurately tracking ionization fronts with the error no larger than one cell. 
Conservation of photons is observed to be at a level of a few percent. The accuracy of the 
results depends weakly on the size of the time step in all time step comparable to the 

typical size of a time step for a cosmological hydrodynamic simulation is sufficient. We also 
apply the algorithm to a realistic density field produced by a cosmological hydrodynamic 
simulation and find that conservation of photons is observed at a level of a few percent. 
These extensive tests indicate that the algorithm is also very accurate. 

The tests performed is based on a first implementation of the basic algorithm outlined 
in the paper. It is clear that there is significant room for further improvement over the 
specific implementation. The current implementation is on a uniform mesh but we think 
that variant approaches (such as with adaptive mesh refinement, for example) based on this 
method are possible and will be explored in the future. 



The work is supported in part by grants NAG5-2759, AST91-08103 and ASC93-18185. 
This implementation of the algorithm will be fine tuned, documented and made available 
in due course. 
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